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ABSTRACT 


Dwarf carbon (dC) stars, main sequence stars showing carbon molecular bands, are enriched by 
mass transfer from a previous asymptotic-giant-branch (AGB) companion, which has since evolved to 
a white dwarf. While previous studies have found radial-velocity variations for large samples of dCs, 
there are still relatively few dC orbital periods in the literature and no dC eclipsing binaries have yet 
been found. Here, we analyze photometric light curves from DR5 of the Zwicky Transient Facility 
for a sample of 944 dC stars. From these light curves, we identify 34 periodically variable dC stars. 
Remarkably, of the periodic dCs, 82% have periods less than two days. We also provide spectroscopic 
follow-up for four of these periodic systems, measuring radial velocity variations in three of them. 
Short-period dCs are almost certainly post-common-envelope binary systems, since the periodicity is 
most likely related to the orbital period, with tidally locked rotation and photometric modulation 
on the dC either from spots or from ellipsoidal variations. We discuss evolutionary scenarios that 
these binaries may have taken to accrete sufficient C-rich material while avoiding truncation of the 
thermally pulsing AGB phase needed to provide such material in the first place. We compare these dCs 
to common-envelope models to show that dC stars probably cannot accrete enough C-rich material 
during the common-envelope phase, suggesting another mechanism like wind-Roche lobe overflow is 
necessary. The periodic dCs in this paper represent a prime sample for spectroscopic follow-up and 
for comparison to future models of wind-Roche lobe overflow mass transfer. 


Keywords: Carbon stars (199), Chemically peculiar stars(226), Close binary stars (254), Common 


envelope evolution (2154), Spectroscopy (1558), Period search (1955) 


1. INTRODUCTION 


Carbon (C) stars are those that show molecular ab- 
sorption bands of C, such as C2, CN and CH, in their 
optical spectra (Secchi 1869). Intrinsic C stars have ex- 
perienced C enrichment via dredge-up of fusion prod- 
ucts from their cores. During the thermally pulsing 
(TP) phase of the asymptotic giant branch (AGB), shell 
He flashes cause strong convection in the shell regions 
bringing C into the atmosphere — the third dredge-up. 
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If the C/O ratio increases above unity, a giant C star 
is formed, since C preferentially binds with oxygen to 
form CO, leaving excess C to form the C molecules of 
C», CN and CH. Thus, it was long thought that all C 
stars were giants on the TP-AGB. 

This made it quite surprising when Dahn et al. (1977) 
found the first main-sequence C star, G77-61. This 
dwarf carbon (dC) star cannot have yet experienced fu- 
sion to create C enhancement, or the third dredge-up, as 
it is a main-sequence star. Dahn et al. (1977) put forth 
a few explanations of this C enhancement on the main- 
sequence, with the favored being that G77-61 was in a 
binary system and had experienced C enhanced mass 
transfer from a previous AGB companion. This for- 
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mer AGB companion would since have become a white 
dwarf (WD) and cooled until no longer detectable along- 
side the dC. Indeed this seemed to have been confirmed 
when Dearborn et al. (1986) found G77-61 to be a radial 
velocity (RV) binary with an orbital period of 245.5 d. 

Today, G77-61 is no longer the only known dC, with 
close to 1000 known in the literature. The majority of 
these dCs come from the Green (2013) and Si et al. 
(2014) samples which were found from all-sky spectro- 
scopic surveys. This has included almost a dozen *smok- 
ing gun" systems in which the WD companion is suffi- 
ciently hot to be visible in the optical spectra as a spec- 
troscopic composite (Heber et al. 1993; Liebert et al. 
1994; Green 2013; Si et al. 2014). These samples have 
shown that the dC stars are actually the most common 
type of C star in the Galaxy. 

Their carbon-enriched atmospheres make dCs the 
most likely progenitors of the carbon-enhanced metal- 
poor, CH, and possibly the Ba II stars, all showing car- 
bon and s-process enhancements (Jorissen et al. 2016; 
De Marco & Izzard 2017). These stars, being more lu- 
minous than dCs, have been studied via RV campaigns, 
which have shown increased binarity compared to nor- 
mal O-rich stars, indicating they have likely also expe- 
rienced mass transfer from an unseen companion (Sper- 
auskas et al. 2016). Barium dwarfs and CH subgiants 
show periods from RV analysis of 1-20 years (Escorza 
et al. 2019). 

Blue straggler stars are another class similar to dCs 
in that they may have experienced mass transfer from a 
previous AGB companion. As discussed by Gosnell et al. 
(2019), blue straggler stars in a cluster color-magnitude 
diagram are more luminous and bluer than the main se- 
quence turnoff. While some are likely formed in mergers 
or collisions, most blue straggler stars are thought, like 
dC stars, to be the result of mass transfer from a gi- 
ant to a main sequence dwarf. Most blue straggler stars 
are found in wide binaries with periods of order 1000 
days, consistent with expectations of mass transfer from 
an AGB star onto a main-sequence companion (Chen & 
Han 2008; Gosnell et al. 2014), which leaves a CO-core 
WD remnant (Paczyński 1971). Those blue straggler 
stars that form after mass transfer from an red giant 
branch star, yields a blue straggler in a shorter binary 
period (of order 100 days; Chen & Han 2008) leaving a 
He-core WD companion. The salient point relevant to 
dC stars is that significant mass is typically gained in 
these encounters. 

While dC stars are now known to be numerous, de- 
tails of their properties remains sparse. This is espe- 
cially true of their orbital properties. Currently, there 
are only six orbital periods for dCs in the literature. 


The first is of the dC prototype G77-61, found to be a 
single-line spectroscopic binary with an orbital period of 
245.5 d (Dearborn et al. 1986). The central source of the 
Necklace Nebula was found to be a binary with a dC, 
which has a photometric period of 1.16 d (Corradi et al. 
2011; Miszalski et al. 2013). The three longest period 
dCs in the literature are those from Harris et al. (2018) 
who found astrometric binaries with periods of 1.23 yr, 
3.21 yr, and 11.35 yr. Margon et al. (2018) found a dC 
with a photometric period of 2.92 d and confirmed this 
as the orbital period with spectroscopic follow-up. 

In their recent work, Whitehouse et al. (2021) con- 
ducted an RV survey of seven dCs with Ho emission, 
finding short orbital periods for all of them (six new 
periods). In addition, they found photometric periods 
with similar lengths as the orbital periods in the range of 
0.2-5.2d. Their light curve modeling suggests that the 
source of variability in their dCs is from stellar rotation 
and spots. As with the new photometrically periodic 
dCs in this paper, these dCs must have experienced a 
common-envelope phase with the former AGB compan- 
ion. 

There have also been large sample few-epoch spec- 
troscopy campaigns of dCs. Whitehouse et al. (2018) 
conducted a few-epoch survey of 28 dCs finding RV vari- 
ability in 21 of them, implying a high binary fraction. 
Roulston et al. (2019) conducted a larger survey of 240 
dCs using few-epoch spectroscopy from the Sloan Digi- 
tal Sky Survey (SDSS; Blanton et al. 2017). They found 
that dCs are consistent with 10096 binarity, with sepa- 
rations of order lau and periods of order lyr. Both 
the Whitehouse et al. (2018) and Roulston et al. (2019) 
surveys lacked enough spectral epochs to fit individual 
orbits and instead relied on statistical analysis to de- 
scribe the dC population as a whole. 

de Kool & Green (1995) modeled the space density 
of dCs, and they predicted the dC period to be bi- 
modal with peaks near LO? LO d and 10?—10? d, con- 
sistent with the known periods listed above. de Kool & 
Green (1995) also found that the production of dCs is 
strongly dependent on metallicity, finding no dCs should 
be formed in systems with initial metallicity greater than 
half of solar (i.e. [Fe/H]> —0.3). This is in agreement 
with metallicity measurements of G77-61, where Plez & 
Cohen (2005) found [Fe/H] = —4.0, making G77-61 one 
of the lowest metallicity stars known. This also has been 
supported by Farihi et al. (2018) who found 30-6096 of 
dCs to be halo objects, which are metal poor. 

Until this year, the known dC periods spanned from 
~ 1d to ^ 41000, likely indicating different formation 
pathways. The longest dC periods that are of order 
10s of years likely experienced only standard Roche-lobe 
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overflow (RLOF) or wind-RLOF (WRLOF). These peri- 
ods are consistent with other types of post-mass-transfer 
systems, such as the blue straggler stars. 

The dCs with periods S 1d would have likely ex- 
perienced common-envelope (CE) evolution (Paczynski 
1976; Ivanova et al. 2013), since the TP-AGB envelope 
expands to several 100s of solar radii. Of interest is how 
CE evolution affects dC formation. Once the CE phase 
has started, the plunge-in of the lower-mass companion 
(in our case, the future dC) could truncate the evolution 
of the AGB by ejection of its envelope. If this happens 
before the TP-AGB phase and the third dredge-up, it is 
likely that the C enhancement needed for dC formation 
will not occur. However, if the CE begins after the AGB 
companion has already become a C-giant, then it may 
be possible for the main-sequence companion to accrete 
enough C-rich material from the CE to become a dC (de- 
pending on the accretion efficiency). If the accretion ef- 
ficiency is not high enough, however, the main-sequence 
companion will not accrete enough material from the 
CE alone, requiring some combination of CE evolution 
with efficient mass transfer before the CE phase that 
is sufficient to transform an O-rich main-sequence star 
into a C-rich dC. 

Significant accretion of mass and angular momentum 
from the AGB companion could result in significant 
spin-up and subsequent activity in some dCs (Green 
et al. 2019). If there are dCs left in very tight orbits 
with the WD remnant, they may show tidally locked 
rotation periods (synchronous rotation), as well as tidal 
distortions causing ellipsoidal variations in photometric 
light curves. A search for periodicity in photometrically 
variable dCs could reveal some systems useful for con- 
straining their evolution. 

Another motivation to study variability in dCs is that 
no dC masses have yet been measured, because there 
are no known eclipsing dC systems. We can estimate 
dC masses from optical or infrared (IR) colors (see Sec- 
tion 6.2), but these estimates have uncharacterized sys- 
tematics, due to differences between normal O-rich stars 
and C stars in the optical and IR regions. 

'This lack of eclipsing dC systems highlights the im- 
portance of photometric surveys to search for the first 
well-characterized eclipsing dC systems. These systems 
could, when combined with RV follow-up, provide the 
first reliable dC mass measurements, and help us un- 
derstand more about the amount, and composition, of 
accreted mass needed to form a dC. 

Margon et al. (2018) conducted a search for periodic 
dCs using the Palomar Transient Factory (PTF; Law 
et al. 2009; Rau et al. 2009), finding just one periodic dC. 
However, they clearly highlighted the potential for large 


photometric surveys to find periodic dCs, particularly 
dCs with short periods that should have experienced 
the strongest phases of CE mass transfer. 

In this paper we report on a unique sample of close 
binary dCs — implicating them as post-common enve- 
lope binaries (PCEBs) and likely pre-CVs — discovered 
from their periodic photometric variability in the Zwicky 
'Iransient Facility. In Section 2 we describe the sample 
of dCs we selected to search and in Section 3 we de- 
tail our process for cleaning and preparing the raw light 
curves. In Section 4 we describe our process for finding 
which dCs have detected periodic signals. In Section 5 
we present spectroscopic follow-up for four of the peri- 
odic dCs in this paper. Finally, in Section 6 we present 
comparisons of these short period dCs to binary popu- 
lation synthesis models to understand how a common- 
envelope phase relates to dC formation. 


2. SAMPLE SELECTION 


'To search for variability in as many dCs as possible, 
we compiled a list of all dCs from the current litera- 
ture. The largest contributor (747 dCs, 7996) is the 
Green (2013) sample of carbon stars from the SDSS. 
We also selected a smaller number of dCs from Si et al. 
(2014), who found 96 new dCs using a label propa- 
gation algorithm from SDSS DRS8, and from Li et al. 
(2018) who selected carbon stars from the Large Sky 
Area Multi-Object Fiber Spectroscopic Telescope sur- 
vey (LAMOST; Cui et al. 2012) using a machine learn- 
ing approach. Our resulting final sample consists of 944 
dCs. 

With our compiled sample, to ensure that any periodic 
candidate was indeed a dwarf carbon star, we used Gaia 
EDR3 parallaxes, proper motions (Gaia Collaboration 
et al. 2021) and distances (Bailer-Jones et al. 2021). We 
required that each periodic C star had Mg > 5mag 
from Gaia EDR3 based either on (1) significant parallax 
w/Wer > 5 (27/34 periodic dCs) or (2) a significant 
proper motion (s/o, > 5) which sets an upper limit on 
the dC distance by limiting its transverse velocity to be 
less than an assumed Galactic escape velocity (Smith 
et al. 2007) of about 600 km/s (7/34 periodic dCs). 


3. LIGHT CURVE PROCESSING 


Using our list of dCs, we cross-matched our sample to 
the Zwicky Transient Facility DR5 (ZTF; Bellm et al. 
2019; Masci et al. 2019; Graham et al. 2019). We re- 
quired a match to be within 2” of our target coordinates 
and each star having > 10 epochs in the available ZTF 
filters. 

From the resulting matches detected within each fil- 
ter, we grouped all sources within the match distance to 
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Table 1. Light Curve Statistics light curve. This variable mean magnitude can cause 
issues with our period search. Therefore, we removed 


these long term trends by fitting out a third-order poly- 


Filter  Nstars < Nepochs > ON < mag ^ < Omag > 


epochs 


ZTF g 833 185 204 19.32 0.11 nomial to the raw light curve. 
ZTE r 867 269 237 18.07 0.05 
ZTF i 554 31 22 17.81 0.05 


4. PERIODIC VARIABILITY 


For each light curve, we searched for periodic signals 


NOTE—Statistics of the light curves in the three ZTF filters. For 
each filter we report the number of stars, the mean number of 
epochs, the standard deviation of the number of epochs, the mean 


magnitude, and the mean magnitude error. 


ensure all epochs for each dC were included. The final 
sample of light curves resulted in 833 dCs with ZTF g 
light curves, 867 dCs with ZTF r light curves, and 554 
dCs with ZTF i light curves. For each light curve, we 
only used epochs for which the ZTF flag catflags —— 
(no ZTF flags), ensuring every epoch is of a high quality. 
We summarize the light curve sample for each filter in 
Table 1. 


We checked for any epochs which appear to be dis- 
crepant by performing an outlier removal on all the 
light curves. We first select from the raw light curve 
the brightest and faintest 596 of epochs. Within these 
brightest and faintest 596, we calculate the median mag- 
nitude of each (i.e. the median of the 5% brightest and 
596 faintest) and the mean error of that same bright- 
est and faintest 596. We then removed any outliers that 
were 2c brighter than the median of the brightest 596, 
and removed those 20 fainter than the median of the 
faintest 596. If this selection dropped the number of 
epochs below 10, we removed that light curve from our 
analysis. This treatment rejects most artifacts without 
removing genuine astrophysical variability. 

We checked the light curves for each dC, in each fil- 
ter, to determine if each dC had detected variability 
by examining how the mean magnitude changed across 
the observed light curve time span. A small number of 
dCs which show no periodic variability in our analysis 
in Section 4 (and a few periodic dCs) show signs of non- 
periodic variability, as well as secular, long-term trends. 
'These non-periodic but variable dCs are of interest and 
may be signs of flaring, variable obscuration, or perhaps 
accretion onto the WD companion. They warrant fur- 
ther investigation, but we do not discuss them further 
in this paper. 

The light curves that show long-term trends of bright- 
ening or dimming on 100s of days timescales cause the 
mean magnitude to vary over the entire time span of the 


down to a minimum period of 0.1 d using the the Lomb- 
Scargle periodogram (LS; Lomb 1976; Scargle 1982). We 
used the Astropy (Astropy Collaboration et al. 2018) 
implementation of the LS algorithm (VanderPlas et al. 
2012; VanderPlas & Ivezić 2015). We selected the high- 
est peak, and if this peak corresponds to an observa- 
tional alias (1d, 29.5d, lyr, etc.) or a harmonic of one 
of these aliases (1/2, 1/3, 1/4, 1/5, 2, 3, 4, 5), we re- 
moved that signal from the light curve and recalculated 
the periodogram until the highest-power frequency was 
not an alias (we counted a frequency not as an alias if it 
was more than 150 frequency bins away from the pure 
alias frequency, i.e. more than 0.005 d-1 away from an 
alias). 

For the highest remaining peak, we calculated the 
false-alarm probability (VanderPlas 2018). We required 
that log (FAP) < —5 in at least one filter for us to select 
a specific dC as a periodic candidate, more conservative 
than e.g., the log(FAP) € —3 used in the recent ZTF 
periodic variable catalog of Chen et al. (2020). 

For the dCs which have light curves selected as peri- 
odic candidates, we checked for any possible harmonic 
confusion in the found period. For each dC, in each 
filter, we plot a power spectrum from the LS analysis. 
'This is used to determine how strong the highest-power 
frequency is compared to the log (FAP) limit and the 
background peaks. Figure 1 shows an example power 
spectrum for an object with a very strong periodic sig- 
nal and shows clear peaks (with 1-d aliasing) above the 
background, and the resulting phased light curve. The 
complete figure set (90 figures) is available in the online 
journal. 

Fig. Set 1. Periodic dC candidate light curves 
and power spectra 

In some cases, the strongest peaks were aliases, typi- 
cally harmonics of 1 month, that overwhelmed the power 
spectrum. For these dCs, we inspected each power spec- 
trum in conjunction with phased light curves. If an- 
other non-alias peak (i.e., with a frequency more than 
0.005 d7! away from an alias) was found in the power 
spectrum meeting our FAP limit, that new peak was se- 
lected as the period for that dC. If no non-alias peaks 
could be found, the dC candidate was removed from our 
sample. 
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RA: 15h19m05.96s DEC: +50d07m02.90 | P = 0.302356d | Amp = 0.106 


Residual 


Frequency [d !] 


Figure 1. Phased light curves and power spectra for all periodic dC candidates. This example light curve and power spectrum 
is for the dC SDSS J151905.964-500702.9. This dC shows a clear and strong signal at 3.3d~* (with 1-d aliasing) that stands out 
above the low noise background in the power spectrum. The red horizontal line represents the peak height needed for a signal 
to meet our log (FAP) € —5 criterion. Grey vertical dashed lines mark the 1-d aliases caused by gaps in data collection, and 
the best period from our analysis is marked by an arrow. The highest significance peak is used to fold the observed light curve, 
yielding the phase-folded light curve in the top panel. Each light curve is plotted twice in phase to clearly show the periodic 
variability. The data are shown as the black scatter points with their respective error bars, and the best fitting model (see 
Section 4) is marked by the red solid line. The residuals are shown below the light curve. The complete figure set (35 figures) 


is available in the online journal. 


Some dCs show strong periodic signals in one filter, 
but do not reach our FAP limit in the other available 
filters. For these dCs, if one filter has a period that 
meets our FAP limit and that period is visible in the 
other filter, we include that second filter even if its FAP 
does not meet our limit. This makes it possible for some 
dCs to have a log(FAP) > —5 in a filter if they have 
log (FAP) € —5 in another filter. 

For all periodic dC candidates selected after inspection 
of their power spectra, we plotted phased light curves 
folded on the highest selected peak period. In addition, 
we plotted 8 different harmonics of that period (1/2, 
1/3, 1/4, 1/5, 2, 3, 4, 5) to check for aliases caused by 
gaps in the observational coverage. Using this plot, we 
calculated model-fit statistics (x?) and selected which 
period harmonic has the best model fit. We used the 
best period to phase the light curve, to which we fit the 
final periodic model. 

Our best-fit models were computed using the auto- 
matic Fourier decomposition (AFD) method, as de- 
tailed in Torrealba et al. (2015). We set an up- 
per limit on the number of Fourier series terms of 
Nmax = 6 to reduce over-fitting. No significant non- 
harmonic terms were included; though one dC, LAM- 
OST J062558.334-023019.4, showed different peaks in 


its power spectrum between the g and r filters with the 
second highest peak in each filter being the highest peak 
in the other. The best AFD model was used to calculate 
the amplitude and epoch of brightest time (to) for each 
light curve. We removed any dC for which the folded 
light curve shows no clear periodic signal or for which 
the amplitude of the variability was less than 0.005 mag. 
For dCs with multiple filters, we use the period from the 
filter with the strongest detection as the selected period 
and force the fits in the other filters to this fixed pe- 
riod. Some filters may not have a clear detection from 
the signal found in another filter, resulting in model fits 
with large errors for that filter. Figure Set 1 contains 
the folded light curves, models, and power spectra for 
all the periodic dC candidates. 

'Table 2 contains the properties for this final periodic 
dC sample. We estimated errors for the best found 
period using a Markov-Chain Monte Carlo (MCMC) 
method. For each dC, in each filter with a detected pe- 
riod, we started 50 MCMC walkers in a Gaussian around 
the detected period. We sampled the walkers for 10,000 
steps each, at each step using the phase dispersion min- 
imization technique (Stellingwerf 1978) to calculate the 
likelihood at each walker position. We used the 1c of 
the marginalized period distribution as the photometric 
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period error for that dC. However this is only a statis- 
tical error, and does not account for the possibility that 
we have selected an alias rather than the true period. 

Our final dC sample contains 34 individual dCs that 
are periodic in at least one ZTF filter. Given the wide 
initial orbits necessary for progenitor dC systems to 
avoid truncation of the TP-AGB phase before enough 
C-rich material can be transferred, it is remarkable that 
19 (5696) of these dCs have periods < 10 (and 28 (8296) 
of these dCs have periods < 2 d), indicating they should 
have experienced a common-envelope (CE) phase. The 
likely origins of the variability in these dCs include spot 
rotation on the dC or tidal distortion of the dC atmo- 
sphere from being in a close orbit with a WD. Since 
many of these dCs have short periods, we assume that 
these systems would have experienced a CE phase and 
have circularized and synchronized (Hurley et al. 2002). 
However, if the light curve variability is from the dC be- 
ing tidally distorted, our detected period would be half 
the orbital period (even with 2x longer true orbital pe- 
riods, these systems should still have experienced a CE 
phase). 

In addition, the 1d aliasing caused by the observing 
window function causes peaks at frequencies of +1d7!. 
These alias peaks can also meet our log(FAP) limit 
(see Figure 1), and while we take the highest signifi- 
cance peak from the filter which produces the best fit- 
ting model (via a x? fit) there is a possibility this is 
the wrong period. This can only be solved by either 
low cadence photometry or confirming the photometric 
period with RV follow-up. For example, the dC SBSS 
1310+561 has +1d-1 aliases the meet our log (FAP) 
limit. In our initial search using ZTF DR3 data, the g 
and r filters had different highest peaks, with the best 
period in the g filter being —5.18d and the best period 
in the r filter being ~0.838d. These are separated ex- 
actly by the 1 d^! aliasing of the window function, with 
the r filter providing a better model fit. However, us- 
ing the newest (and larger) ZTF DR5 data set results in 
both the g and r filters having the same highest peak, at 
5.1878 + 0.0012d. Whitehouse et al. (2021) confirmed 
this as the period with their RV observations of this dC. 

The recent work by Whitehouse et al. (2021) included 
modeling of light curves for their sample of periodic dCs. 
They examined whether spot rotation, tidal distortion 
or irradiation by a hot WD companion could be the 
source of the photometric variability in their dC sam- 
ple. They found that for periods near and longer than 


1d, both tidal distortion and irradiation are reduced to 
levels that would not be detectable in the light curves. 
While tidal distortion could be detectable for our short- 
est period dCs in this paper (0.1 « P « 0.2d), the 
predicted amplitudes at these periods from Whitehouse 
et al. (2021) are larger than those we find in our sam- 
ple (as was the case for their sample). The irradiation 
modeling of Whitehouse et al. (2021) (which assumes 
a WD temperature from 30,000 K to 20,000 K) predicts 
amplitudes large enough to be detected. However, the 
majority of our dC amplitudes are smaller than pre- 
dicted by the models, suggesting irradiation is not the 
source of variability for most of our dCs. Additionally, 
the majority of the dCs in our sample are mid- and late- 
type dCs (see Roulston et al. 2019) which do not have 
a visible WD in their optical spectra. This sets a limit 
that for these types, that WD must be cooler than about 
10,000 K, reducing the irradiation effects below our de- 
tection limits. Only the composite dC+WD systems 
which have a hot WD (like SDSS J151905.964-500702.9) 
may have detectable irradiation effects. This leaves spot 
rotation as the most likely source of variability in our pe- 
riodic dC sample. However, the origin of the variability 
in these dCs is not truly confirmed without comparison 
to spectroscopic RV follow-up. 

Finally, as the ZTF survey continues to accure more 
data we expect to find more photometrically variable 
dCs from the current sample of known dCs. However, 
even given a favorable inclination (say i > 85?) the ZTF 
errors are too large for the detection of an eclipse of a 
cool WD in these systems. Using our estimated dC radii 
and luminosities and assuming a WD companion with a 
standard mass of 0.6Mọ and temperature of 7000K (we 
expect to see the WD component in the optical spec- 
tra if it is any hotter than this) we would only expect 
an average primary eclipse depth of 0.005 mag. This is 
below our detection threshold with ZTF, with our dCs 
having median errors of 0.019 mag, compared to their 
median amplitude of 0.059 mag. The Vera Rubin Ob- 
servatory's LSST survey (Ivezić et al. 2019) is expected 
to have errors of approximately 0.005 mag for a point 
source with r — 19.0, which may allow for the detec- 
tion of dC eclipses. However, the majority of known 
dCs reside outside the LSST footprint, with 1796 below 
the declination cut of of ó = +2° (Ivezié et al. 2019). 
Detection and characterization of the first eclipsing dC 
system will likely require dedicated observations with 
high cadence and low photometric noise. 
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5. SPECTROSCOPIC FOLLOW-UP 


To constrain the origins of the photometric variability 
we have begun spectroscopic follow-up of the periodic 
dCs discovered here. We report spectroscopic follow-up 
for four of these dCs: SDSS J151905.964-500702.9, SDSS 
J123045.53--410943.8, LAMOST J062558.334-023019.4 
(referenced further on as J1519, J1230, J0625 respec- 
tively) and SBSS 1310+561. 


5.1. Spectroscopic Set-Up 


The dCs J1519 and J1230 were observed with the Bi- 
nospec spectrograph on the MMT telescope (Fabricant 
et al. 2019). For all observations, we used the 0.85” slit 
with the 600 1 mm `! grating centered on 7250 A, giving 
coverage from 6000 Á to 8000A covering Ha and the 
CN bands. The reduced spectra have a dispersion of 
0.61 A pix 1 with R¥3590. All Binospec data were re- 
duced using the standard Binospec reduction pipeline! 
(Kansky et al. 2019). 

The dC J0625 was observed with the Magellan Echel- 
lette (MagE; Marshall et al. 2008) spectrograph on the 
Magellan Baade Telescope. All observations used the 
0.85” slit and were reduced using the MagE reduc- 
tion pipeline? (Chilingarian 2020). The reduced spectra 
cover from about 3200 Á to 10000 A with R~4500. 

Observations for SBSS 1310+561 were acquired at 
the 1.5m Fred Lawrence Whipple Observatory (FLWO) 
telescope with the FAST spectrograph (Fabricant et al. 
1998) using the 600 1 mm^! grating and the 1.5"slit, 
which provides wavelength coverage from 6000À to 
8000 A at 1.5 À spectral resolution. 


5.2. SDSS J151905.96+500702.9 


One of the more interesting dCs in our periodic sam- 
ple, with photometric periodicity detected with highest 
significance, is SDSS J151905.964-500702.9 (also known 
as CBS 311; we use J1519 in the rest of this paper), a 
dC+DA spectroscopic composite binary. J1519 was dis- 
covered by Liebert et al. (1994) and has been studied 
on numerous occasions (Farihi et al. 2010; Green 2013; 
Whitehouse et al. 2018; Ashley et al. 2019; Roulston 
et al. 2019; Green et al. 2019). However, this is the first 
reporting of its periodic variability. 

J1519 (r — 17.3 mag) has four epochs of optical spec- 
train the SDSS, with the most recent spectrum shown in 
the top panel of Figure 2. The spectrum of J1519 shows 
a dC with a hot DA WD companion, as well as Ho emis- 
sion. Whitehouse et al. (2018) and Roulston et al. (2019) 


l https:/ /bitbucket.org/chil.sai/binospec/wiki/Home 
? https:/ /bitbucket.org/chil_sai/mage-pipeline/src/master/ 


found RV variability using few-epoch spectroscopy with 
ARVmax of 46.8 + 15.8 km s-! and 44 + 20 km sl, 
respectively. Farihi et al. (2010) conducted a study of 
WD red dwarf systems, including J1519, using the Hub- 
ble Space Telescope. They found 01519 to be unresolved, 
placing the constraint on its separation of < 10 au. 


5.2.1. J1519 WD Model Fits 


Since J1519 is a spectroscopic dC+DA composite, we 
can fit WD model atmospheres to the WD component 
to fit Tog and log(g) using the SDSS spectra. Bédard 
et al. (2020) fit WD models and found fit values of 
312304210 K and 7.97+0.05 respectively. 

Farihi et al. (2010) found that spectroscopically fit 
WD parameters are often biased due to a cool com- 
panion. To update the fits of Bédard et al. (2020), we 
performed our own model atmosphere fits to the DA 
component of J1519 using the synthetic WD model at- 
mospheres of Levenhagen et al. (2017). We first fit the 
late type dC (dCM) template of Roulston et al. (2020) 
to the SDSS spectrum of J1519 by finding the best-fit 
velocity, shifting the template, and then scaling it to the 
flux near Ha. We then removed the dC spectrum from 
the total spectrum, leaving just the WD component. We 
then fit the visible Balmer lines from HO and blue-ward 
to the entire grid of WD model spectra. We interpo- 
lated the grid of WD model spectra to include half-steps 
in the model space. Our best-fitting model parameters 
for Teg and log(g) were 31000+500 K and 7.85+0.05, re- 
spectively, and can be seen in Figure 2. The black line is 
the single SDSS spectrum with the highest S/N shifted 
to the rest-frame, and the blue line is the best fit WD 
model spectrum. We did not use Ho for the WD fit as 
the dC component contributes most to the spectrum in 
emission. In addition, we did not use the H9 line, as 
only half of the line is visible in the SDSS spectrum. 

'The WD temperature of our fit is in good agreement 
with that of Bédard et al. (2020). However, our log(g) is 
0.12 dex lower, resulting in both our WD mass and cool- 
ing age being lower than those in Bédard et al. (2020). 
For the purposes of this paper, we adopt our fit values 
of log(g) and Tez. The WD properties we use can be 
found in Table 3, with the mass, radius, and cooling age 
coming from the models of Fontaine et al. (2001). 


5.2.2. J1519 Radial Velocities 


Although RV variability has been detected in J1519, 
there are no published RV orbital fits for this system. 
Based on our photometric analysis, we found a period 
of 0.302356 + 0.000021 d (— 7 hr) for J1519. Therefore, 
we conducted a spectroscopic monitoring of J1519 using 
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Table 3. J1519 WD Properties 


Parameter Value Error Source 


Tore [K] 31000 500 (1) 

log g [dex] 7.85 0.05 (1) 

M [Mo] 0.57 0.02 (2) 

R [Re] 0.015 0.001 (2) 

Tcoo1 [Myr] 7.7 0.2 (2) 
NoTrE—Best fit model parameters 
for the DA component of SDSS 
J151905.964-500702.9. Each parameter 
lists the source used: (1) this paper (2) 


from evolutionary models of Fontaine 
et al. (2001) 


the MMT spectroscopic setup as was described in Sec- 
tion 5.1. On the nights of 2020 August 19 and 20, we 
observed a sequence of 21x 200 s exposures, on the night 
of 2020 August 22 we observed 27x 2005s exposures, and 
on the nights 2021 April 21 and 23 we observed 24x230 s 
exposures. The exposures on each night were then co- 
added in threes, resulting in seven final epochs on the 
first two nights, nine epochs on the third night, and eight 
on each of the last two nights for a combined total of 39 
epochs (with about 600s total exposure each), with an 
average S/N = 5 for all epochs in the continuum region 
near Ha. 

Since the full spectrum includes both stellar compo- 
nents, we measured the RV from the Ha emission line, 
presumed to come from the dC atmosphere alone. First, 
for each epoch, we re-scaled the late-type (dCM) dC 
template of Roulston et al. (2020) to the flux in our 
MMT spectrum in the region of 6300-6500 À. We then 
used this as the model for the dC continuum level of that 
epoch, which was used to calculate the Ha emission line 
center, equivalent width and associated errors. The RV 
measurements have an average error of approximately 
5 km s^! , and the equivalent width measurements have 
an average error of approximately 0.18 À. 

Figure 2 shows the measured RV (middle) and Ha 
equivalent widths (bottom) for J1519. 'To fit the RV 
curve, we used the rvfit program which uses a simu- 
lated adaptive annealing procedure, the details of which 
can be found in Iglesias-Marzoa et al. (2015). We left 
all parameters free to be fit, with the solution quickly 
converging to a circular orbit. We therefore refit the RV 
curve leaving all parameters free again except for the 
eccentricity, which we fix to e = 0.0. The resulting best- 
fit model can be seen in Figure 2 (red curve) and the fit 
parameters can be found in Table 4. 

The best fitting orbital period from the RVs 
(0.327526 + 0.000012 d) is longer than the best photo- 


metric period by 0.025170 + 0.000024d (about 36 min- 
utes). Fixing the period in the RV fitting procedure 
to that of the photometric period results in a poorer 
model fit, with the longer period model being a bet- 
ter fit at the 3.2c level. The best-fit semi-amplitude of 
K> = 33.314 km s^! (ARVmax = 2K = 66.6 km s^!) 
is in agreement with the RV variations found by White- 
house et al. (2018) and Roulston et al. (2019), as their 
random epochs likely did not catch the true RV ampli- 
tude. However, the low measured semi-amplitude sug- 
gests an extremely low inclination of this system, with 
i & 10? if we take our estimated dC mass of 0.41 Mo 
from Section 6.2. 

One possible explanation for a longer orbital period 
than photometric period is that J1519 was spun up by 
the accretion that it experienced, and has not yet syn- 
chronized the rotation and orbital periods in the approx- 
imate 8 Myr since mass transfer stopped (assuming the 
mass transfer ceased at the same time the WD formed). 
Green et al. (2019) analyzed Chandra observations of 
J1519 (as well as five other Ha emission dCs) and found 
it to show X-ray emission consistent with having a short 
rotation period, which would lend support to the accre- 
tion spin-up scenario. Deeper photometric imaging and 
RV follow-up, particularly of the WD component, could 
even better characterize this system. It is clear, however, 
that this dC has both photometric and RV variability 
on a <0.33-d timescale, indicating it most likely has a 
short orbital period and formed through a CE event. 


5.3. SDSS J123045.538+410943.8 


Another interesting dC is SDSS J123045.534-410943.8 
(J1230), whose SDSS spectrum is shown in Figure 3. 
J1230 shows the Cə and CN lines typical of late type 
dC stars, but also shows strong absorption lines of K 
and a strong CaH band near 6800 A. Additionally, this 
dC shows strong emission lines of Ha, H8, Hy, Hó, and 
Ca H and K. Unlike J1519, there is no visible WD com- 
ponent in the spectrum of J1230. 

We observed 01230 (r = 16.6 mag) with the Binospec 
spectrograph on the 6.5m MMT using the setup de- 
scribed in Section 5.1. We took 6266s spectra on the 
nights of 2021 February 3, 5, 7, 8, 9, 11 and 15. This 
resulted in a total of 42 spectra, with an average S/N 
of 13 in the continuum region near Ha. We measured 
the line center and equivalent widths of the Ha line, as 
well as for the two K lines visible in our MMT spectra of 
J1230. The emission and absorption lines have the same 
velocities, indicating the are coming from the same re- 
gion. We average these velocities to measure the RV for 
each of the 42 epochs. 
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Figure 2. TOP: SDSS spectrum for the dC+DA SDSS J151905.96+500702.9, in black. The hot WD is visible, as are the 
carbon bands of C2 and CN. The spectrum also shows Ha emission. The best fitting model atmosphere to the DA component 
of the dC-- WD composite J1519 is shown in blue, with the inset to the right showing a zoomed-in, stacked view of the Balmer 
lines used in the model fit (the black lines are the normalized flux of the WD component, and the blue lines are the best fitting 
model atmosphere, dashed grey lines means that that Balmer line was not used in the fit due to poor quality). The resulting 
dC component — simply the observed spectrum with the WD model subtracted — is shown in red. MIDDLE: RV as measured 
from the Ha line for J1519, phased on the time of periastron passage from the RV fit, and the fit RV period of 0.327526 d. The 
red solid line represents the best fitting model. BOTTOM: Equivalent widths measured from the Ho line, phased on the fit RV 
period of 0.327526 d. The blue curve is the best-fit model to the data of a single sinusoid. The y-axis has been inverted so that 


smaller equivalent width values (more emission) are up. 


In the same method as J1519, we used the rvfit pro- 
gram to fit the RV curve of J1230. For this dC, we left 
all parameters free for fitting, with the resulting best 
period fit matching that of the photometric light curve 
(P = 0.882519 + 0.000020 d). We therefore fix the pe- 
riod to the photometric period and the eccentricity to 
0.0, and refit the RVs. The resulting best fit can be 
found in Table 4 and the phased RV curve in Figure 3. 
As with J1519, we find the parameter errors using a 
MCMC method centered around the best fit parame- 
ters. 

The best fit gives a circular orbit with a semi- 
amplitude of K> = 123.0 + 0.7 km s71. If we use our 
estimated dC mass from Section 6.2 (0.25 Mọ) and an 
assumed WD mass of 0.6 Mo the implied inclination of 
this system is around 7 = 56°. 

'The presence of multiple emission lines of H and Ca 
suggest that the photometric variability of J1230 is com- 
ing from re-processing of the WD flux on the surface of 
the dC. Even though the WD companion to J1230 is not 
hot enough to be seen in the optical spectrum, it may 
be warm enough to still heat the surface of the dC. If 


this is true, we could expect the dC to be at maximum 
brightness when the WD-facing side is pointed toward 
us maximally, i.e. when the dC is moving transversely 
on the sky, between the ascending and descending nodes. 
Comparing the light curve of J1230 to the RVs however 
shows this is not the case, as if it were, we would expect 
the RV to be moving to the descending node after the 
photometric maximum, which the RV curve for J1230 
is 0.33 out of phase with. This suggests that the pho- 
tometric variability is not coming from re-processing, 
but rather from spot rotation on an active dC (with the 
emission lines indicating chromospheric activity). We 
do note that the uncertainties on our epoch of maxi- 
mum brightness and period may cause us to incorrectly 
predict the phase when our spectroscopy was collected 
in 2021 February by up to 0.03 cycles. 


5.4. SBSS 1310+561 


We observed SBSS 13104-561 (r — 14.1 mag) using the 
FAST spectrograph on the 1.5m telescope at FLWO us- 
ing the setup described in Section 5.1. We took 3x300s 
spectra during the nights of 2021 February 10 and 11, 
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Figure 3. TOP: SDSS spectrum for the dC SDSS J123045.53+410943.8. The typical carbon bands of Cz and CN for late 
type dCs are labeled. Additionally, clear and strong emission of Ho, H8, Hy, Hó, and Ca H and K are visible. Middle: RV as 
measured from the Ho and K I lines for SDSS J123045.534-410943.8, folded at the photometric period of 0.882519 d. The red 
dashed line represents the best fitting model. Bottom: Equivalent widths measured from the Ha line, phased on the photometric 
period. The blue curve is the best fit single sinusoid model to the data, while the grey dotted line is the average equivalent 
width value. The y-axis has been inverted so that smaller equivalent width values (more emission) are up. We do not detect 
any significant variation in phase for the equivalent width, to a limit of «1.50 À. 


and 6x300s spectra on the night of 2021 February 12 
for a total of 12 spectra with an average S/N of 32 in 
the continuum region near Ha. 

Unlike with our MMT and Magellan observations, 
because of observing time constraints on our awarded 
FAST time, we chose to obtain these spectra close to 
the quadrature phases based on the ZTF photometry 
(P = 5.1878 + 0.0012d). We assumed that the pho- 
tometric period corresponds to the orbital period, and 
used to from the light curve to calculate the expected 
times that SBSS 13104-561 should be at the quadrature 
phases (¢ = 0.25 and $ = 0.75). Our actual observations 
were taken at phases 0 = 0.27 -E0.02 and 9 = 0.47-E0.01. 

From these spectra, we measured the RV at 0 = 0.27 
to be V, = —79.7 + 9.5 km s^! and the RV at 0 = 0.47 
to be V. = —19.3 + 5.5 km s^!. Taking the difference 
in these two velocities for this system (ARV = 60 + 
11 km 8 1) can place a lower limit on the semi-amplitude 
of Kə > 30 + 6 km el Using our estimated mass 
from Section 6.2 (0.46M) and an assumed WD mass 
of 0.6 Mo, this constrains the inclination to i > 25° (if 
i = 60°, then we would expect K = 109 km s^). Since 
the phase difference between our two epochs is quite 


small, this ARV suggests that SBSS 1310+561 is in a 
tight orbit and is very likely a PCEB. 


5.5. LAMOST J062558.33--023019.4 


The dC LAMOST J062558.33--023019.4 (hereafter 
J0625, r — 13.9 mag) was observed on the nights of 2021 
January 11 and 12 using the Magellan MagE instrument 
setup described in Section 5.1. Each night, we observed 
15x300s exposures. The final reduced spectra consists 
of 30 epochs with an average S/N of 22 each in the con- 
tinuum region near Ha. 

Using the Ha emission line, we measured the RV of 
J0625 for each epoch. We found no evidence for RV 
variability, nor any variability in Ha equivalent width. 
We found the RV to vary with only a standard deviation 
of 3.9 km s 1, and with a maximum ARV = 12.1 + 3.2 
km al. In addition, cross-correlation of the spectra 
across epochs resulted in no significant measured RV 
variations. Using our estimated mass from Section 6.2 
(0.84 Moe) and an assumed WD mass of 0.6 Mç, this 
places a constraint on the inclination of i € 7? (if i — 
60°, then we would expect K> = 44 km s^!). This may 
suggest that the photometric variability is not related 
to the orbital period in this system since such a low 
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Table 4. Radial Velocity Fits 


Parameter J1519 J1230 
P [d] 0.327526 + 0.000012 0.882519” 
T, [MJD] 59080.2085 + 0.0053 — 59265.07955 + 0.00059 
e 0.0" 0.0" 
w [deg] 90.0" 90.0° 
y [km s -1] -1.7 + 2.8 -2.9 + 0.5 
Kə [km s^ !] 33.3 + 14 123.0 + 0.7 
a3 sini [Re] 0.2153 + 0.0090 2.15 + 0.01 
f(mı, m2) [Mo] ` 0.00125 + 0.00016 0.170 + 0.003 
xà 2.6 0.97 
Nobs 39 42 
'Time span [d] 247.2 12.15 
rmsa [km s TI 11.6 2.5 


? Parameter fixed during fitting. 


NOTE—Fit parameters from the radial velocity follow-up. The value 
for each parameter is given as the median of the marginalized dis- 
tribution of the MCMC samples. The errors for each parameter are 
the Lo values from the marginalized distribution of the MCMC sam- 
ples. Additionally, derived values for the orbital separation (asin i) 
and mass function (f (mi, m2)) are given. 


inclination (and low semi-amplitude) is unlikely if the 
photometric period of 7.6080 + 0.0014 d represents the 
orbital period. Hence, this system adds weight to the 
evidence that the photometric variability in dCs may 
often be due to spot rotation. 


6. COMMON ENVELOPE CONNECTION 


For the progenitor of the dC companion to become a 
C giant, it must enter the third-dredge up phase (Iben 
1974). AGB stars have 8 degenerate CO core with a 
double-shell (moving outward from the core) of helium 
and hydrogen. As the hydrogen shell (which produces 
most of the energy) continues to fuse H into He, the 
helium shell surrounding the core continues to grow. 
Eventually, the helium shell experiences runaway fusion, 
driving expansion of the envelope material above. This 
He-shell “flash” and expansion means the star is now 
in the thermally pulsing AGB (TP-AGB) phase. He- 
lium shell fusion causes the inter-shell region to become 
strongly convective, dredging helium fusion products to 
the surface, i.e., the third dredge-up. As the expansion 
continues, the pressure in the helium shell will drop, 
eventually stopping its energy production. The layers 
contract again with hydrogen shell fusion resuming, and 
the cycle repeats. 

Each successive thermal pulse becomes stronger, 
reaching deeper into the intershell zone, and the stel- 
lar radius increases (Iben & Renzini 1983). As helium 
shell fusion products are brought to the surface, it is 


possible that the envelope carbon abundance increases 
until C/O > 1. Since C preferentially binds with O, C2 
and CN bands only appear when C/O> 1, forming a C 
giant star. 

AGB stars going through the TP-AGB phase can 
reach radii of 800 Re (3.7 au) as they experience succes- 
sively stronger thermal pulses (Marigo et al. 2017). As- 
suming an AGB mass of 2.5 Mo, AGB radius of 800 Ro, 
and a dC mass of 0.4 Mo, this system would experience 
the beginning of a common-envelope (CE) with an ini- 
tial period of z 4.2 yr (if the dC mass is 1.0 Mo instead, 
then Pæ 3.8yr). Therefore, dCs with initial periods 
~ Ayr (1500 d) or less will very likely have experienced 
a CE phase, corresponding to the shorter-period peak 
modeled by de Kool & Green (1995). The dCs in this 
paper with P« 1d are most certainly the result of a CE 
spiral-in. Of the six dC periods in the current literature, 
two of them have P« 3d, so have likely experienced a 
CE. It seems then that many dCs may have experienced 
8 CE phase. 

Dell'Agli et al. (2021) recently studied the extreme 
AGB stars (those AGB stars which have extremely red 
mid-IR colors, e.g. Gruendl et al. 2008) and showed 
that the excess dust and outflow densities of these stars 
may be explained by envelope stripping in a CE event. 
'Their models suggest that these extreme AGB stars are 
actually post-common-envelope binaries (PCEBs) with 
orbital periods of order 1 d, matching the periods for dCs 
in our sample. Dell’Agli et al. (2021) also found that 
the CE in their models starts after the rapid growth 
of the AGB radius, once the C/O ratio increases past 
unity, which corresponds well with the requirements for 
producing the short-period dCs we find. This makes 
these extreme AGB stars potential progenitors systems 
of the dCs that are in the CE phase currently. 

However, is mass accretion during a CE phase the 
most likely mass transfer mechanism to form dCs? We 
can address this question by looking at our periodic dC 
sample in the context of models that simulate expected 
binary populations. 


6.1. Binary Population Synthesis Models 


We used the binary population synthesis (BPS) mod- 
els of Toonen & Nelemans (2013) to see if the observed 
population of dCs can be reproduced by theory. The 
full details of the BPS models can be found in Toonen 
& Nelemans (2013) and are briefly described here. 

These BPS models were created using the SeBa 
(Portegies Zwart & Verbunt 1996; Nelemans et al. 2001; 
Toonen et al. 2012; Toonen & Nelemans 2013) popula- 
tion synthesis code. This code generates an initial pop- 
ulation of binaries and simulates their evolution, taking 
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into account processes such as stellar winds, magnetic 
breaking, mass transfer, common-envelope, and angular 
momentum loss. The initial stellar population is gener- 
ated from the classical BPS distributions found in Too- 
nen & Nelemans (2013) via a Monte Carlo method. The 
resulting binaries are then convolved with a Galactic 
model including a star formation history that depends 
on time and location in the Milky Way based on Boissier 
& Prantzos (1999) so that the simulated binaries can be 
compared to our observed sample. 

For the synthetic populations used here, the common- 
envelope phase is modeled on the basis of the energy 
budget i.e. the classical a-formalism of Tutukov & Yun- 
gelson (1979). We discuss the results of two different 
models here that account for two different CE efficien- 
cies: model aa and oo? which have «A of 2 and 0.25, 
respectively. The parameter A is the structure param- 
eter of the envelope to calculate the envelope binding 
energy (Paczynski 1976; Webbink 1984; de Kool et al. 
1987; Livio & Soker 1988; de Kool 1990; Xu & Li 2010). 
The o parameter describes the efficiency with which or- 
bital energy is consumed to unbind the CE. A smaller 
value of o implies less efficient usage of orbital energy, 
and therefore a stronger shrinkage of the orbital period 
during the CE-phase. We do not consider the orbital 
angular momentum method of Nelemans et al. (2000), 
as this model does not reproduce the observed charac- 
teristics of the general PCEB (WD /main sequence) pop- 
ulation (Toonen & Nelemans 2013). 

Furthermore, the BPS models here allow for accre- 
tion during the common-envelope phase. The accretion 
rate is limited by the thermal timescale of the accretor 
times a factor that is dependent on the stellar radius and 
the corresponding Roche lobe (Portegies Zwart & Ver- 
bunt 1996; Toonen et al. 2012) following Kippenhahn 
& Meyer-Hofmeister (1977); Neo et al. (1977); Packet 
& De Greve (1979); Pols & Marinus (1994). The total 
accreted mass is then given by the integral of the ac- 
cretion rate times the timescale of the CE event, which 
here is taken to be 100 yr. This timescale is consistent 
with hydrodynamical simulations (Ricker & Taam 2008; 
Ivanova et al. 2013) and observations of hot subdwarf 
binaries (Igoshev et al. 2020), although cataclysmic vari- 
ables may suggest a longer CE timescale, up to 10* yr 
(Michaely & Perets 2019; Igoshev et al. 2020). 


6.2. BPS Comparison to Observed dC Sample 


We use the resulting model population for a direct 
comparison to our observed sample of short-period dCs, 
assuming the photometric period is the current orbital 
period. To do this, we estimated dC masses based on 
their infrared absolute magnitude Mx in the K band. 


Comparisons of M dwarf spectra (Ivanov et al. 2004) 
to C star spectra (Tanaka et al. 2007) reveal them to 
be much more similar in the infrared than in the op- 
tical region. We used K. band (2.159 um) magnitudes 
from the Two Micron All-Sky Survey (2MASS; Skrut- 
skie et al. 2006). Six of our periodic dCs do not have K, 
band magnitudes. For these, we first fit Gaia absolute 
G band (Mg) to the dCs that do have K, band magni- 
tudes. This fit was then used to convert the Gaia Mg 
into Mx, for the dCs lacking K, band magnitudes. We 
then fit Mx, for our dCs to stellar masses using data 
from Kraus & Hillenbrand (2007). This fitting also pro- 
vides us with bolometric luminosities for the dCs in our 
sample. Comparing our bolometric luminosities to those 
provided in Green et al. (2019) (who used a spectral en- 
ergy distribution method fitting 0.35 — 12.5 um) for the 
four dCs that overlap, we find our luminosities agree 
within 3%, indicating our dC mass estimates should be 
reliable. 

Our mass estimates can be found in Table 5. We find 
that none of our dCs are fit with masses > 1M or 
< 0.2Mo, in agreement with the range for which de- 
tectable C2, CN, and CH bands are expected. We note 
that some of the lowest mass dCs may have been brown 
dwarfs or even planets before they accreted significant 
C-enriched material from their former AGB companion. 

Using the mass-radius relationship for main-sequence 
stars of Eker et al. (2018), we estimate the radius for 
these periodic dCs as well, which are included in Table 5. 
Using these estimated radii we calculate the Roche-lobe 
filling factor (RLFF), using the equation of Eggleton 
(1983) to find the Roche radii. Six out of 34 of our 
periodic dCs may be experiencing RLOF back onto the 
WD (all have a RLFF > 1 in Table 5). However, we 
caution that physical parameters are derived from O- 
rich main-sequence models, which may not accurately 
represent all dCs. For example, (1) we do not know 
the mass of the unseen WD companion and assume it 
is 0.6 Mo (2) we assume these mass-radius and Mx- 
mass relations hold for dCs, as they do for normal O- 
rich stars (3) dCs are thought to be of a lower metallicity 
population and studies have found that low metallicity 
M dwarfs may have smaller radii (Kesseli et al. 2019) and 
(4) since dCs may have increased activity and magnetic 
fields due to their mass accretion, their radii may be 
inflated (Kesseli et al. 2018). We see no obvious evidence 
of flickering or accretion outbursts in any of our ZTF 
light curves that might indicate current RLOF back onto 
the WD. 
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Table 5. Periodic dC Parallaxes, Distances, and Estimated Physical Parameters 


R.A. (J2016.0) Decl. (2016.0) ^ «^ ot. a" ei  BP-RP' MG Mk Mac logig(Lpg/Lo) Raq  RLFF 
[mas] ^ [mas] [pc] [pc] [mag] [mag] [mag] [Mo] [erg s] [Ro] 
00h47m06.76s  +00d07m48.80s 0.68 0.18 1340 196 178 785 Am 0.67  —0.88 0.59 0.06 
Olh3lm19.05s  +37d20m25.30s 1.12 — 0.12 944 88 159 763 479 uer  —0.88 0.59 0.27 
02h35m30.65s +02d25m18.58s 1.078 0.08 590 30 161 793 5.08 0.60 1.05 0.52 0.22 
02h54m14.24s +26d21m54.19s 3.294 0.082 — 301 8 146 725 452 0.74 —0.72 0.07 1.410 
O4hl6m05.lls  +50d28m28.52s 2.946 0.015 335 2 1.19 604 4.00 0.87 0.40 0.83 0.12 
05h02m40.82s +40d23m23.59s 1.287 0.086 840 62 161 7.58 465 0.71  —0.79 063 0.13 
0Gh25m58.34s +02d30m19.43s 2410 0.024 — 409 4 1.09 593 3.90 0.90  —0.34 0.86 0.11 
OTh44m4T.60s +51d38m31.76s 2.178 0.050 457 11 164 792 5.09 0.60 1.05 0.52 0.23 
O8hlim57.14s  +14d35m33.00s 1.596 0.089 612 14 0.69 672 439 0.77 0.64 0.70 0.45 
O9hl4m58.08s +21d56m39.65s 3.594 0.050 275 4 1.82 843 523 0.56 —1.14 0.49 0.25 
09h33m24.58s ^ -O0d31m44.07s 5.726 0.036 173 1 1.63 791 514 0.59  —1.08 0.81 0.27 
09h40m26.288 +36d25m48.81s 1.55 — 0.21 765 — 92 1.77 8.93 5.61 048 —1.35 0.41 0.17 
2h02m46.01s  -t54dl9m29.24s 1.08 0.19 1103 170 1.98 8.87 532 0.55  -1.18 047 0.26 
2h08m53.35s -00d08m47.99s 0.78 — 037 243 372 1.35 695 465° 0.71 —0.80 063 0.70 
2h10m06.99s +58d43m18.34s 1.134 0.064 873 44 141 746 4.79 0.67 —O.88 0.59 1.04 
2h23m57.62s +55d01m51.43s 1.911 0.079 521 21 1.79 9.03 5.86 043 —1.48 0.36 0.50 
2h30m45.52s +41d09m43.45s 5.736 0.056 173 1 244 10.38 682 0.25 —1.96 0.22 0.20 
3h03m59.18s  +05d09m38.62s 1.44 0.10 722 58 1.53 7.96 516 0.58  —140 0.50 0.20 
3hl2m42.27s +55d55m54.84s ^ 9.54 0.023 106 1 191 910 571 046  -—140 0.39 ^ 0.08 
13h31m23.61s  --48d26m24.37s 1.10 0.15 959 136 1.92 8.92 5.62 ous  -135 0.40 0.75 
4h09m53.08s  -06dllm4l.71s 2.502 0.079 393 11 1.32 645 411 0.84  -0.47 0.79 0.87 
4h15m15.24s  +51041028.016 0.77 0.36 4420 799 119 6.02 447° 0.75  —0.69 0.68 — 0.87 
5hlim44.58s +38d59m10.46s 2.05 — 0.11 487 29 184 8.90 588 042  -1.49 0.35 0.50 
5hl5m42.72s +52d01m45.47s 1.291 0.065 775 37 1.61 8.00 530 0.55  -1.18 0.47 0.60 
5hl9m05.93s +50d07m03.14s 2274 0064 437 14 0.77 9.08 5.95 O41  -153 0.34 0.52 
5h24m34.12s +44d49m55.84s 0.96 0.19 1236 211 2.00 8.89 — 5.80 — 044  -1.45 0.37 0.62 
5h25m04.49s +32d25m10.90s 0.10 0.40 3280 598 1.55 742 491° 0.64  —0.95 0.56 — 1.22 
5h30m59.26s +45d12m00.338s ^ 2.215 0.044 444 8 181 845 539 0.53  -122 0.45 0.05 
5h35m32.92s +01d10m16.22s 1.18 0.39 1179 234 2.15 9.09 480 0.67  -—0.88 0.59 1.07 
Gh37m18.63s  +27d40m26.63s ^ 2.497 0.077 397 12 242 9.54 5.922 041  -151 0.35 0.21 
Gh59m02.30s +25d05m49.00s 0.88 0.24 1293 220 2.07 887 5.79 — 044 — —144 0.37 0.57 
9h23m55.93s +44d58m32.20s 1.238 0.048 802 28 1.36 6.90 479 067  —0.88 0.59 121 
22h08m10.015 +25d17m30.17s 4.121 0.053 241 3 1.56 7.75 4.79 uer — —0.88 0.59 ` 0.60 
23h4lm30.74s +15d19m43.20s 0.27 0.12 2685 482 1.04 5.07 399 0.87 0.40 0.83 1.60 


? From Gaia Collaboration et al. (2021) 
b From Bailer-Jones et al. (2021) 


^ Mk interpolated from MG 


NorE—Distances and magnitudes for the periodic dCs in our sample. We use the Gaia distances, colors, and magnitudes, as well as the 2MASS 
absolute K magnitudes to estimate masses and bolometric luminosities for our dCs. For the solar bolometric luminosity, we adopt the value 
log, Lo = 33.58. We also calculate the Roche-lobe filling factor (RLFF) under the assumption of a 0.6 Me WD companion. We calculate the 
mass errors to be of order 0.05 Mo, the logy, (tel! Lo) errors to be of order 0.1, and the radius errors to be of order 0.05 Ro. However, we caution 


that physical parameters are derived from O-rich main-sequence models, which may not accurately represent all dCs. 


To compare the BPS models directly to our observed 
dC sample, we applied a series of cuts and selection ef- 
fects to the models, as follows: (1) P < 100d (2) r mag 
< 19.5 (3) Mac < 1Mo (4) 1.0 <Mzams < 4 Mo (5) 
the initial primary must be a TP-AGB star at the on- 
set of the CE phase. Here, Mac is the current mass of 
the main sequence companion in the BPS models, and 
Mzams is the initial mass of the primary at the begin- 
ning of the models (which will become the AGB donor). 


We show the resulting BPS models in Figure 4 and 
Figure 5 — models oo? and aa, respectively. In both 
figures, the BPS models are shown as the colored 2D 
histogram in mass and period (note that the histogram 
color scale is logarithmic and its range is different for 
each plot), and the periodic dC stars from this paper 
represented as the red scatter points (with KDE con- 
tours). The dashed black line represents the RLOF 
boundary, with systems occupying the region to the left 
(shaded in grey) filling their Roche lobes, under the as- 
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Figure 4. Binary population syntheses results for model «o2, which is the model with lower common-envelope efficiency 
(aà = 0.25) in which less orbital energy goes to unbinding the CE. The background heat map shows the number of systems 
(colored on a log scale) in model aa2 going through a CE phase while the primary is in the TP-AGB phase. The plotted masses 
are for the secondary (i.e. the dC) and periods are for the currently observed system. Each panel uses a different magnetic 
braking formalism: (a) magnetic braking from Rappaport et al. (1983) (b) magnetic braking from Ivanova & Taam (2003) (c) 
magnetic braking from Knigge et al. (2011). The red plus scatter markers are our observed short period dC sample without 
Ha emission, while the purple triangles are dCs which show Ha emission. The red contours are a KDE contour for the entire 
periodic dC sample. The dashed black curve represents the boundary for current dC systems to be experiencing RLOF, where 
systems to the left in the grey shaded region may be experiencing RLOF. The bottom panel shows a histogram of the estimated 
mass accreted by the secondary during the CE phase. While this model (aa2) reproduces the period distribution better than 
model aa (see Figure 5), the estimated accreted mass is even lower than that of aa because of the lower CE efficiency. Since 
20.03 Mo of mass transfer is likely required to create a dC, a CE is not likely to be the primary mechanism for accretion to 
form a dC. 


dard magnetic braking of Rappaport et al. (1983). From 
Figure 5, it is seen that this model is unable to repro- 


sumption of a 0.6 Mọ WD companion. Both figures also 
show a histogram of the estimated mass accreted during 


the CE phase (assumed to last 100 yr). 

Figure 4 shows model «a2 (aA = 0.25) and includes 
three different magnetic braking prescriptions. Panel (a) 
uses the magnetic braking of Rappaport et al. (1983), 
panel (b) that of Ivanova & Taam (2003) and panel (c) 
that of Knigge et al. (2011). Again, the color scale is 
logarithmic and its range is different for each sub-figure. 

Model aa2, however, does not reproduce the mass dis- 
tribution of our dCs very well, generating low mass sys- 
tems than observed (still under the assumption that our 
physical parameters derived from O-rich main-sequence 
models apply to dCs). While it may be that model oo? 
does not produce low mass dCs, we have not considered 
our sample selection effects in this comparison. Our ob- 
served sample is likely biased toward lower mass dCs as 
(1) they have stronger C9 and CN bands, and (2) their 
variability is fractionally larger and so easier to detect. 

Model oa (Figure 5) uses a higher CE efficiency (aA = 
2) similar to classical BPS studies and includes the stan- 


duce the short periods of our observed dC sample. This 
is in agreement with the conclusions based on the SDSS 
PCEBs (WD4-MS systems; Zorotovic et al. 2010; Too- 
nen & Nelemans 2013; Camacho et al. 2014), where Too- 
nen & Nelemans (2013) found that standard efficiency 
(aA = 2) CE was also unable to reproduce the observed 
periods, as it generated too many long-period PCEBs. 

A crucial shortfall is that the estimated mass accreted 
for all models is too small to convert a main-sequence 
star into a dC (see Section 7 for a discussion). Mis- 
zalski et al. (2013) suggest that to shift the secondary 
envelope from approximately solar (C/O); — 1/3 to the 
observed (C/O); > 1 requires AM» = 0.03-0.35 Mo for 
M» = 1.0-0.4 Mc. The predicted mass accretion in our 
BPS models is lower than this by 2—3 orders of magni- 
tude. Together with the strong mismatch between the 
modeled and observed dC period-mass distributions, it 
seems clear that there must be further mass accretion 
outside the brief CE phase. 
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Figure 5. Binary population synthesis results for model 
aa, which is the model with higher CE efficiency (aA = 2). 
The background heat map shows the number of systems (col- 
ored on a log scale) in model aa going through a common- 
envelope phase while the primary is on the TP-AGB phase. 
The plotted masses are for the secondary (i.e. the dC) and 
periods are for the currently observed system, including stan- 
dard magnetic braking from Rappaport et al. (1983). The 
red plus scatter markers are our observed short period dC 
sample without Ha emission, while the purple triangles are 
dCs which show Ha emission. The red contours are a KDE 
contour for the entire observed dC sample. The dashed black 
curve represents the boundary for current dC systems to be 
experiencing RLOF, where systems to the left in the grey 
shaded region may be experiencing RLOF. The bottom panel 
shows a histogram of the estimated mass accreted onto the 
secondary during the CE phase. This model (ao) does not 
reproduce our observed dC period distribution, as it does not 
produce periods below 1d, as well as not being able to ac- 
crete enough mass (70.03 Me) we expect necessary to form 
dCs. 


Qualitatively it is also possible to argue that accre- 
tion during CE evolution is rarely significant for non- 
degenerate companions. The common envelope itself 
typically possesses much higher specific entropy than 
the surface of the accretor, with the consequence that 
matter accreted by the companion star reaches pressure 
equilibrium at the surface of that star with much higher 
temperature, and vastly lower density, than the accre- 


tors initial surface layer. A temperature inversion or 
roughly isothermal layer is expected to bridge this en- 
tropy jump with the result that, over the duration of 
the CE evolution (which is much shorter than the ther- 
mal time scale of the accretor), the accretor is thermally 
isolated from the common envelope, while the common 
envelope itself becomes increasingly tenuous. 

A solution to the under-prediction of accreted mass 
may be that more mass accretion may take place be- 
fore the CE phase, but during the TP-AGB phase, by 
wind accretion during wind-RLOF (WRLOF; Mohamed 
& Podsiadlowski 2007). WRLOF is a mass transfer 
state that lies between standard wind mass transfer and 
standard RLOF, where the wind of the primary star 
is focused toward the secondary star. This results in 
increased mass transfer efficiency as compared to the 
standard Bondi-Hoyle-Lyttleton case (Hoyle & Lyttle- 
ton 1939; Bondi & Hoyle 1944). 

In the WRLOF regime, the primary is technically not 
filling its Roche lobe. However, low-velocity wind mat- 
ter is funneled through the Roche lobe to the companion, 
allowing for mass transfer to take place in binaries with 
wider orbits than the classical RLOF case. WRLOF 
would boost dC formation since, if the initial orbital 
separation is too small, the expanding AGB atmosphere 
can cause a CE before the third dredge-up can turn the 
AGB into a C star. 

A variety of simulations (Abate et al. 2013; Saladino 
et al. 2018, 2019; Saladino & Pols 2019) have shown 
that WRLOF in binaries with AGB primaries can have 
mass-transfer efficiencies of 40-5096. For an average 
AGB wind mass loss rate of 1077-1074 Maer") (Hófner 
2015), a main sequence companion could accrete enough 
material (~ 0.35M5) in only 105—108 yr, within the time 
an AGB star is expected to stay a C star (109yr; Marigo 
et al. 2017). WRLOF has also been shown to efficiently 
tighten the orbit (Saladino et al. 2018; Chen et al. 2018) 
so that more systems could be driven towards orbits 
with the periods we find. 

Indeed, it appears the WRLOF may be the dominant 
mass transfer mechanism for many chemically peculiar 
stars. Abate et al. (2013) showed that simulations of 
AGB binaries with WRLOF were better able to re- 
produce the observed formation rates of the CEMP-s 
stars. Saladino et al. (2019) and Saladino & Pols (2019) 
also performed simulations finding AGB binaries with 
WRLOF were consistent with the observed properties of 
the CEMP-s, CH, and Ba stars. This further strength- 
ens the connection between dCs, WRLOF, and the other 
chemically peculiar stars. 

However, detailed modeling of the WRLOF in the spe- 
cific case of the TP-AGB phase with a C star is needed to 
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understand how the larger stellar radii of AGB-C stars 
and the increased dust formation often found in their 
winds affect the WRLOF mass transfer efficiencies. The 
periodic dCs in this paper represent a prime sample that 
is ready for spectroscopic follow-up and for comparison 
to future models of WRLOF mass transfer. 


7. SUMMARY 


We searched ZTF light curves of a sample of 944 dCs 
for periodic signals. We found 34 dCs with signs of 
significant photometric variability, with 8296 having P« 
2d. The most likely origins of this periodicity is either 
from spot rotation or surface heating of the dC from 
the close WD companion. Even if the detected ZTF 
periods arise from ellipsoidal variations and represent 
half the orbital period, such short periods are surprising 
for dC stars, which require significant accretion from 
a TP-AGB C giant to turn them into the C-enriched 
dwarfs we see today. 

Spectroscopic follow-up is needed to determine the 
source of the variability in each of these periodic dCs, 
especially to confirm for the case of spot rotation that 
the system is circularized and tidally locked so that the 
rotation period can be assumed to equal the orbital pe- 
riod (ie. that our reported photometric periods cor- 
respond to both rotational and orbital periods). The 
periodic dCs in this paper provide a rich new sample to 
target for spectroscopic follow-up, as well as to study 
dC formation and properties. We have confirmed the 
photometric period as the orbital period for one dC for 
which we have obtained spectroscopic follow-up. In two 
other dCs we have confirmed that they must have short 
orbital periods from their RVs (not confirming the pho- 
tometric period as the orbital period however). In all 
three cases, these short (P« 1d) orbits indicate these 
dCs have indeed experienced a CE phase. 

'These short periods indicate that at least some dCs 
will experience a CE at some point in their formation, 
with P« 1d dCs having experienced experiencing sub- 
stantial plunge-in. We used binary population synthesis 
models to show that the observed sample of dCs is not 
well-reconstructed by mass transfer during the common- 
envelope phase alone, since the dCs in our sample re- 
quire at least 0.03 Mọ of mass accretion but our mod- 
els predict 2-3 orders of magnitude less transfer dur- 
ing the CE phase, suggesting mass accretion before the 
CE phase. However, some systems such as cataclysmic 
variables indicate CE timescales an order or two longer 
(Michaely & Perets 2019; Igoshev et al. 2020) than our 
assumed 100 yr (based on Ricker & Taam 2008; Igo- 
shev et al. 2020), which may substantially increase the 


amount of accreted material to the point that the CE 
alone could provide enough mass to form a dC. 

Hydrodynamical simulations of CE evolution typically 
find that accretion onto a non-degenerate companion is 
not common, because of the entropy barrier between the 
companion and the surrounding material (e.g. Ivanova 
et al. 2013, for a review). In fact, even in the case of neu- 
tron star companions, which can accrete more efficiently 
due to neutrino cooling, accretion is limited to < 0.1 Mo 
(MacLeod & Ramirez-Ruiz 2015). Further modeling of 
the CE phase involving C-AGBs may provide further 
insight. 

dC systems that begin as very wide binaries would 
experience stable mass transfer and a widening of the 
orbit. Systems that initially are close would begin a 
CE phase either during the red giant branch or during 
the AGB before the TP-AGB and, without experienc- 
ing the third dredge-up during the TP-AGB, would not 
produce dCs. Therefore, it seems that the most likely 
mass transfer mechanism to form dCs is WRLOF. 

Further modeling of WRLOF in binaries with a TP- 
AGB star are needed to fully test this formation pathway 
of dCs. Additionally, further work is needed to under- 
stand the relationship between initial dC metallicity and 
mass to constrain the amount (and composition) of ma- 
terial that needs to be accreted to form a dC. This would 
be an important step in constraining the mass-transfer 
efficiency in the WRLOF case. 


Facilities: FLWO:1.5m (FAST), Magellan:Baade 
(MagE), MMT (Binospec) 


Software: Astropy (Astropy Collaboration et al. 
2013, 2018), Corner (Foreman-Mackey 2016), Matplotlib 
(Hunter 2007), Numpy (Harris et al. 2020), Scipy (Vir- 
tanen et al. 2020), Scikit-learn (Pedregosa et al. 2011), 
TOPCAT (Taylor 2005) 
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